Fusariose: QTLs Description - Part1

1/ Introduction

This file aims to describe the QTLs found concerning fusariose resistance.
It concerns the evaluation of every blups computed before.
Let’s upload the file containing all the LOD scores. It contains about 1.6M lines, so it is a bit long

#Watch out, to reproduct analysis, you have to update the path.
data=read.table("~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/bilan_simple_marker.gz" , header=T , sep="," )
#data=data[ order(data$LG , data$Distance) , ]
# random lines of the dataset to develop quicker
#data=data[ sample( seq(1,nrow(data)) , 30000 )  , ]
#data=data[ order(data$LG , data$Distance) , ]

Let’s upload the genotyping matrix:

# Genotype matrix
geno = read.table(file="~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/fichier_genotypage_QTL.csv", sep = ";" , header = F, na.strings = "-")
geno= as.matrix(geno)
colnames(geno)= geno[1,]
geno= as.data.frame(geno[-1 , ])
# Remove useless individuals
#geno <- geno [ geno[,1] %in% names(data),]

Let’s upload the phenotyping matrix:

#Watch out, to reproduct analysis, you have to update the path.
pheno=read.table("~/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/PHENOTYPE/phenotypage_all_fusa_blup.csv" , header=T , sep=";" )
#numeric
pheno[,-1]=apply(pheno[,-1],2,as.numeric)
#check no redondancy in genotype name
table(pheno$geno)[table(pheno$geno)>1]
## named integer(0)
# put geno name as rowname
rownames(pheno)=pheno$geno
pheno=pheno[,-1]
# delete columns with only NA
which(apply( pheno , 2 , function(x) all(is.na(x)) )==TRUE)
## named integer(0)
pheno=pheno[ , ! apply( pheno , 2 , function(x) all(is.na(x)) ) ]

Charge some libraries that will be useful

library(RColorBrewer)
library(rmdformats)
library(xtable)
library(plotly)
library(FactoMineR)
library(knitr)

Let’s create a color vector for the study

my_colors = brewer.pal(8, "Set2") 
my_colors = colorRampPalette(my_colors)(20)

How many variables do we have?

my_carac=levels(data$variable)
nlevels(data$variable)
## [1] 115

What is the LOD threshold to declare a QTL significant? I declare it to 3.73

my_QTL_thres=3.73

Let’s do some object that contains group of variable:

# reps
REP1=my_carac[grep("rep1", my_carac )]
REP2=my_carac[grep("rep2", my_carac )]
BLUP=my_carac[grep("blup", my_carac )]

# experiments
CAP13=my_carac[grep("CAP13", my_carac )]
GRI11=my_carac[grep("GRI11", my_carac )]
GRI13=my_carac[grep("GRI13", my_carac )]
GRI15=my_carac[grep("GRI15", my_carac )]
LEC14=my_carac[grep("LEC14", my_carac )]
BEQ11=my_carac[grep("BEQ11", my_carac )]
BAR14=my_carac[grep("BAR14", my_carac )]

# Type de phéno NEPI, NEPIL, PEPIL, NOT, DON
PEPIL=my_carac[grep("PEPIL", my_carac )]
PEPI=my_carac[grep("PEPI", my_carac)]
PEPI=setdiff(PEPI, PEPIL)
NEPIL=my_carac[grep("NEPIL", my_carac )]
NOT=my_carac[grep("NOT", my_carac )]
DON=my_carac[grep("DON", my_carac )]

# AUDPC et derniere notation
AUDPC=my_carac[grep("AUDPC", my_carac )]
LAST=c("CAP13_PEPI550_rep1","CAP13_NEPIL550_rep1","CAP13_PEPIL550_rep1","CAP13_PEPI550_rep2","CAP13_NEPIL550_rep2","CAP13_PEPIL550_rep2","CAP13_PEPI550_blup","CAP13_NEPIL550_blup","CAP13_PEPIL550_blup","GRI11_NOT550_rep1","GRI11_NEPIL550_rep1","GRI11_NOT550_rep2","GRI11_NEPIL550_rep2","GRI13_PEPI550_rep1","GRI13_NEPIL550_rep1","GRI13_PEPIL550_rep1","GRI13_NOT550_rep1","GRI13_PEPI550_rep2","GRI13_NEPIL550_rep2","GRI13_PEPIL550_rep2","GRI13_NOT550_rep2","GRI13_PEPI550_blup","GRI13_NEPIL550_blup","GRI13_PEPIL550_blup","GRI13_NOT550_blup","LEC14_PEPI550_rep1","LEC14_NEPIL550_rep1","LEC14_PEPIL550_rep1","LEC14_PEPI550_rep2","LEC14_NEPIL550_rep2","LEC14_PEPIL550_rep2","LEC14_PEPI550_blup","LEC14_NEPIL550_blup","LEC14_PEPIL550_blup","BEQ11_PEPI450","BEQ11_NOT500","BEQ11_NEPIL500","BAR14_PEPI600_rep1","BAR14_NEPIL600_rep1","BAR14_PEPIL600_rep1","BAR14_PEPI600_rep2","BAR14_NEPIL600_rep2","BAR14_PEPIL600_rep2","BAR14_PEPI600_blup","BAR14_NEPIL600_blup","BAR14_PEPIL600_blup")

# agronomical variable
EPI=my_carac[grep("_EPI" , my_carac )]
FLO=my_carac[grep("FLO" , my_carac )]
HEI=my_carac[grep("HEI" , my_carac )]
VER=my_carac[grep("VER" , my_carac )]

# Resistance set for publication = BLUps when avail. DO we have all the needed traits?
a=c(DON, NOT, PEPI, NEPIL, PEPIL)
b=intersect(a,BLUP)
c=intersect(BEQ11, a)
PUBLI=c(b,c,"LEC14_DON")

Let’s check what is available for which experiment?

# start table
bilan=data.frame( matrix(0,5,7) )
colnames(bilan)=c("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11","BAR14")
rownames(bilan)=c("DON?", "PEPI?", "NOT?", "NEPIL?", "PEPIL?")
all=list(CAP13, GRI11, GRI13, GRI15, LEC14, BEQ11, BAR14)
num1=0
num2=0
for(i in all){
  num1=num1+1
  num2=0
  for(j in list(DON, PEPI, NOT, NEPIL, PEPIL)){
   num2=num2+1
   a=intersect(i, j)
   bilan[ num2 , num1]=length(a)
  }}
kable(xtable(bilan))
CAP13 GRI11 GRI13 GRI15 LEC14 BEQ11 BAR14
DON? 1 0 1 0 0 0 0
PEPI? 4 0 4 3 4 3 5
NOT? 0 6 4 3 0 5 1
NEPIL? 4 6 4 3 4 5 6
PEPIL? 4 0 4 3 4 0 5

We know that:
- DON is available for CAP13, GRI13 and BAR14 only. - PEPI is not available for GRI11 - NOT is not available for CAP13. - PEPIL is not available for GRI11 and BEQ11. So please check it is indeed what we have in this QTL file!

2/ Useful functions

A set of function allowing to characterize significant QTLs for a set of variable:

#Same function but with interactive graphics!
plot_a_interactiveQTL=function(data, select_chromo, my_variable, LOD_threshold, my_ylim){
  data=data
  select_chromo="4B"
  my_variable=HEI
  
  p=data %>%
    filter(LG==select_chromo & variable %in% my_variable) %>%
    group_by(variable, Distance) %>%
    select(Distance, variable, LOD, R2, group_physique, Posi_physique) %>%
    summarise(
      LOD=max(LOD),
      R2=max(R2, na.rm=T),
      Posi_physique=mean(Posi_physique, na.rm=T),
      group_physique=table(group_physique)[1]
    ) %>%
    mutate(my_text=paste("<br>", "Pos phy: ", group_physique , " | ", Posi_physique , sep="")) %>%
    ggplot(aes(x=Distance, y=LOD, color=variable)) +
      geom_line()
 
# interactive
ggplotly(p)
}

#A function to plot a QTL and visualize the LOD scores along the chromosome. Works for several variable on a specific chromosome.
plot_a_QTL=function(data, select_chromo, variable, LOD_threshold, my_ylim){
  
  # order data
  data=data[ order(data$LG , data$Distance) , ]

  # plot the first variable
  a=data[data$variable==variable[1] & data$LG==select_chromo & !is.na(data$LOD) , ]
  plot(a$LOD ~ a$Distance , type="l" , ylim=c(0,my_ylim) , ylab="LOD score" , xlab="position (cM)" , col=my_colors[1] , lwd=1.6 )
  
  # for every other variable, add a LOD score
  num=1
  for( i in c(2:length(variable))){
    num=num+1
    a=data[data$variable==variable[num] & data$LG==select_chromo & !is.na(data$LOD) , ]
    points(a$LOD ~ a$Distance , type="l" , ylim=c(0,my_ylim) , ylab="LOD score" , xaxt="n" , col=my_colors[num] , lwd=1.6 )
    }
  
  # add lod threshold
  abline(h=my_QTL_thres , col="grey")

  # add legend
  legend("topright" , horiz=F , col=my_colors[c(1:length(variable))] , legend=variable , bty="n" , lty=1 , pt.cex=2, lwd=1.6)
  legend("topleft" , horiz=F , col="white" , legend=select_chromo , bty="n" , lty=1 , pt.cex=0, lwd=1.6)

  }

#A function to summarize all the significant QTL of a list of variables. Gives the max LOD score, the IC, the position etc..
summary_table_for_QTL=function(data, selected_variable, LOD_threshold, size_IC){

  # order data
  data=data[ order(data$LG , data$Distance) , ]
  my_colors=brewer.pal(5,"Paired")
  
  # create empty summary table
  bil=data.frame(matrix(0,1000,10))
  colnames(bil)=c("pop","carac","chromo","LOD_max","position","marker","IC","R2","a","diff")
  
  # loop to study every variable & chromosome
  num_line_bilan=0
    num_chromo=0

    for(chrom in levels(data$LG)){
        num_variable=0
        num_col=0
        num_chromo=num_chromo+1

        for(var in selected_variable){
            num_col=num_col+1
            current_data=data[ data$variable==var & data$LG==chrom  , ]
            current_data=current_data[!is.na(current_data$LOD) , ]
            signif_data=current_data[ current_data$LOD > LOD_threshold  , ]
            
            #Si j'ai des marqueurs significatifs
            if(nrow(signif_data)>0){
                
                #Je récupère les infos de ce QTL
                LOD_max=max(signif_data$LOD)
                mark_max=signif_data$marqueur[signif_data$LOD==LOD_max]
                pos_max=signif_data$Distance[signif_data$LOD==LOD_max]
                r2_max=signif_data$R2[signif_data$LOD==LOD_max]
                a_max=signif_data$a[signif_data$LOD==LOD_max]
          my_diff=ifelse( signif_data$moy.A[signif_data$LOD==LOD_max][1] > signif_data$moy.B[signif_data$LOD==LOD_max][1] , "Dic2>Silur" , "Silur>Dic2" )
            
                #Détermination IC --> je bloque la zone max a 30 cM
                in_IC=current_data[ current_data$LOD > (LOD_max-size_IC) & current_data$Distance > (pos_max-30) & current_data$Distance < (pos_max+30)  , ]
                IC_min=min(round(in_IC$Distance,2))
                IC_max=max(round(in_IC$Distance,2))
                
                #Je remplie le tableau bilan
                num_line_bilan=num_line_bilan+1
                bil[num_line_bilan , 1]="DS"
                bil[num_line_bilan , 2]=var
                bil[num_line_bilan , 3]=chrom
                bil[num_line_bilan , 4]=LOD_max
                bil[num_line_bilan , 5]=pos_max[1]
                bil[num_line_bilan , 6]=as.character(mark_max)[1]
                bil[num_line_bilan , 7]=paste(IC_min, IC_max, sep="-")
                bil[num_line_bilan , 8]=r2_max[1]
                bil[num_line_bilan , 9]=a_max[1]
                bil[num_line_bilan , 10]=my_diff
                
                #close 3 loops
                }}}
    
# Clean and print the Table
bil=bil[bil$pop!=0 , ]
bil=bil[ order(bil$pop , bil$chromo , bil$carac) , ]
return(bil)
#close function
}

A set of function allowing to study confidence intervals:

#A function to represent confidence interval of QTLs on the genetic map.
show_IC_QTL_on_map=function(data, selected_variable, selected_chromosome, LOD_threshold, size_IC){

  # user can pick every variable
  if(selected_variable[1]=="all"){selected_variable=levels(data$variable)}
  
  # order data
  data=data[ order(data$LG , data$Distance) , ]
  my_colors=brewer.pal(12,"Set3")
  my_colors = colorRampPalette(my_colors)(30)
  var_with_QTL=c()
  
  # loop to study every variable & chromosome
    num_chromo=0

    # prepare background of chart 
    nchro=length(selected_chromosome)
    xchro=c(seq(5,5*nchro,5) )
  par(mar=c(0,5,2,1))
  plot(1,1,col="transparent",bty="n", xaxt="n", yaxt="l" , xlab="", ylab="position (cM)" , xlim=c(3,max(xchro)+11) , ylim=c(320,-5) )
  num=0
  for(chrom in selected_chromosome){
    num=num+1
    A=data[data$variable==selected_variable[1] & data$LG==chrom, ]
    points(rep(xchro[num],nrow(A)) , A$Distance , pch=20 , cex=0.8 , col="grey")
    }

    # IC calculation and plot
  for(chrom in selected_chromosome){

        num_variable=0
        num_chromo=num_chromo+1
        num_col=0

        for(var in selected_variable){
          num_col=num_col+1
            current_data=data[ data$variable==var & data$LG==chrom  , ]
            current_data=current_data[!is.na(current_data$LOD) , ]
            signif_data=current_data[current_data$LOD>LOD_threshold  , ]
            
            #Si j'ai des marqueurs significatifs
            if(nrow(signif_data)>0){
              
              # I add it to the list of variables with QTL
              if(!var%in%var_with_QTL){var_with_QTL=c(var_with_QTL,var)}
                
                # Je récupère les infos de ce QTL
                LOD_max=max(signif_data$LOD)
                mark_max=signif_data$marqueur[signif_data$LOD==LOD_max]
                pos_max=signif_data$Distance[signif_data$LOD==LOD_max]

                # Détermination IC --> je bloque la zone max a 30 cM
                in_IC=current_data[ current_data$LOD > (LOD_max-size_IC) & current_data$Distance > (pos_max-30) & current_data$Distance < (pos_max+30)  , ]
                IC_min=min(in_IC$Distance)
                IC_max=max(in_IC$Distance)
                
                # J'ajoute le trait corresponant a ma variable
                num_variable=num_variable+0.2
                lines(c(xchro[num_chromo]+num_variable,xchro[num_chromo]+num_variable) , c(IC_min,IC_max) , col=my_colors[which(var_with_QTL%in%var)], lwd=6)
                num_variable=abs(num_variable)

        #Close 2 loops and a if
                }}}
    
#Ajout légende et nom de chromosome?
text(xchro , rep(-10,3) , selected_chromosome , col="orange")
#legend("bottomright", 14 , 240 , horiz=F , col=my_colors[1:length(var_with_QTL)] , legend=var_with_QTL , bty="n" , lty=1 , lwd=6 )

#Close function
}

#A function that gives IC features: size, # of markers, # of genes...
find_IC_infos=function(min_IC, max_IC, chrom, to_remove){

  # Get the part of the genetic map concerned by this IC:
  IC_map=data[ which(data$LG==chrom & data$Distance>=min_IC & data$Distance<=max_IC & data$LG==data$group_physique & data$variable==levels(data$variable)[1]), c(1:5)]

  # remove markers if needed.
  if(length(to_remove)>0){
    IC_map=IC_map[-which(IC_map$marqueur%in%to_remove) , ]
  }
  
  # Calculate statistics

  # genet
  deb_IC_gen=min(IC_map$Distance)
  end_IC_gen=max(IC_map$Distance)
  size_IC_gen=end_IC_gen - deb_IC_gen
  size_IC_gen=round(size_IC_gen , 2)
  nb_of_mark=nrow(IC_map)

  # physic
  deb_IC_phy=min(IC_map$Posi_physique)
  end_IC_phy=max(IC_map$Posi_physique)
  size_IC_phy=end_IC_phy - deb_IC_phy
  nb_of_gene=nrow(POS[which(POS$group==chrom & (POS$position*1000000)>=deb_IC_phy & (POS$position*1000000)<=end_IC_phy) , ])

  # return a vec with all info
  vec=c(chrom, min_IC, max_IC, size_IC_gen, nb_of_mark, deb_IC_phy, end_IC_phy, size_IC_phy, nb_of_gene)
  return(vec)
  }

#A function to show all the genes present in a IC. The relatioship between physical and genetic position is also available
show_genes_of_IC=function(min_IC, max_IC, chrom, my_lwd, my_cex, to_remove){
  
  # get data for the IC
  IC_map=data[ which(data$LG==chrom & data$Distance>=min_IC & data$Distance<=max_IC & data$group==data$group_physique & data$variable==levels(data$variable)[1] & data$LG==data$group_physique ) , c(1:5)]

  # remove markers if needed.
  if(length(to_remove)>0){
    IC_map=IC_map[-which(IC_map$marqueur%in%to_remove) , ]
  }

  # physical informations
  #deb_IC_phy=min(IC_map$Posi_physique)
  #end_IC_phy=max(IC_map$Posi_physique)
  my_genes_pos=POS[which(POS$chromo_BW==chrom) , ]

  #plot
  plot(IC_map$Distance ~ as.numeric(IC_map$Posi_physique/1000000),
     ylim=c(min_IC, max_IC),
     xlab="physical pos. of markers (Mb)", ylab="genetic position (cM)" , 
     main=paste("chromome ",chrom ," -- IC : ",min_IC," - ",max_IC," cM",sep="" ),
     col=rgb(0.4,0.3,0.6,0.4) , pch=20 , cex=my_cex
     )
  abline(v=my_genes_pos[,3]/1000000 , col="grey", lwd=my_lwd)

  }

What happens on a specific position

# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_a_QTL=function(marker, trait, ...){

  # prepare data: merge expression and alleles:
  nucl= geno[,c(1,which(colnames(geno)==marker)) ]
  my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
  colnames(AA)=c("Genotype", "Phenotype", "Allele")
  AA$Allele[which(AA$Genotype=="dic2")]="A"
  AA$Allele[which(AA$Genotype=="silur")]="B"
  AA=na.omit(AA)
  AA$Genotype
  
  # Highlight Dic2 and Silur
  AA$type=factor(ifelse(AA$Genotype=="dic2" | AA$Genotype=="silur","Highlighted",AA$Allele))
  
  # Make plot
  p=ggplot(AA, aes(x=Allele, y=Phenotype, text=Genotype, fill=type, color=type )) +
    geom_boxplot(alpha=0.5) + 
    geom_jitter() + 
    theme(legend.position="none") +
    ggtitle(trait)

  # interactive
  ggplotly(p)
  }

What about epistasic relationships

#A function to calculate interaction between 2 QTLs for a trait
# Function to calculate interaction effect between QTLs:
analyse_inter=function(marker1, marker2, trait){

  # prepare data: merge expression and alleles:
  nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
  my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
    AA=na.omit(AA)

    # Prepare
    my_mark=paste(AA[,3],AA[,4],sep="-")
    means <- round(tapply(AA[,2],my_mark,mean,na.rm=T) ,2)
    
    # Calculate interaction
    model=lm(AA[,2] ~ AA[,3] * AA[,4])
    sum=summary(model)
    tot_r2=round(sum$r.squared,3)
    inter_pval=round(sum$coefficients[4,4],10)
    
    # Complete summary files
    res=c(trait, marker1, marker2, means, tot_r2, inter_pval)
    return(res)
}

#A function to boxplot 2 QTLs together
# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_two_QTL=function(marker1, marker2, trait, ...){

  # prepare data: merge expression and alleles:
  nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
  my_pheno= data.frame(ind=rownames(pheno), pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
    AA=na.omit(AA)

    # Prepare
    my_mark=paste(AA[,3],AA[,4],sep="-")
    means <- round(tapply(AA[,2],my_mark,mean,na.rm=T) ,2)

  #plot
    par( mar=c(5,5,2,2))
  boxplot(AA[,2] ~ my_mark , 
      medlwd=0, cex.axis=0.6, cex.col="grey", las=1, main="" ,
      ...,
      ylab=trait,  xlab="Bi-locus genotype" , col=my_colors[c(3,6,8,3)] , xaxt="n", boxwex=0.4) 
  my_labels=c( expression(paste('Dic2'["1B"],"-",'Dic2'["5A"])) , 
               expression(paste('Dic2'["1B"],"-",'Silur'["5A"])) ,  
               expression(paste('Silur'["1B"],"-",'Dic2'["5A"])) , 
               expression(paste('Silur'["1B"],"-",'Silur'["5A"])) )
  axis(labels=my_labels , at=c(1,2,3,4) , side=T)
}

# A function to represent a QTL. It makes a boxplot with 2 categories: alleles A and B
boxplot_two_QTL_interactive=function(marker1, marker2, trait, ...){

  # example to develop if needed
  #marker1="Cluster_16778|Contig1|original@1840"
  #marker2="Cluster_9940|Contig1|complementarySeq@104" 
  #trait="BAR14_PEPIL300_blup"

  # prepare data: merge expression and alleles + add dic2 and Silur
  nucl= geno[ , c(1, which(colnames(geno)%in%c(marker1,marker2)) ) ]
  qq=data.frame(c("dic2","silur"), c("A","B"), c("A","B"))
  colnames(qq)=colnames(nucl)
  nucl=rbind(nucl,qq)
  my_pheno= data.frame(ind=rownames(pheno), pheno=pheno[ , trait])
  AA=merge(my_pheno, nucl, by.x=1, by.y=1, all=T)
    AA=na.omit(AA)

    # Prepare
    AA$my_mark=as.factor(paste(AA[,3],AA[,4],sep="-"))
    means <- round(tapply(AA[,2],AA$my_mark,mean,na.rm=T) ,2)
  AA=AA[order(AA$my_mark),]
  AA$text=paste("ind: ",AA$ind,"<br>","value: ",round(AA$pheno,2),sep="")

  # First we calculate a new x axis with the jiiter added to it. Jitter level must be relative to the number of points per level.
  my_prop=summary(AA$my_mark) / nrow(AA)
  new_x=c()
  for(i in c(1:nlevels(AA$my_mark))){
    myjitter<-jitter( rep(i, summary(AA$my_mark)[i]), amount=my_prop[i]/2)
    new_x=c(new_x, myjitter)
    }

# Plot
plot_ly(AA) %>%
    add_trace( x=as.numeric(AA$my_mark), y = ~pheno, color = ~my_mark, type = "box" ) %>%
    add_trace( x=new_x, y = ~AA$pheno,  
        type = "scatter", 
        mode="markers",  
        marker = list(size = 10, color = ifelse(AA$ind=="dic2","blue",ifelse(AA$ind=="silur","red",rgb(0.1, 0.1, 0.1, 0.3))), line = list(color = rgb(0.1, 0.1, 0.1, 0.8), width = 2)),
        text=~text,
        hoverinfo="text"
    )%>%
    layout(title = "",
           hovermode="closest",
         showlegend=F,
         xaxis = list(title = "" , tickfont = list( color="orange", size=26 ),  tickmode="array", tickvals=c(1,2,3,4) , ticktext=levels(AA$my_mark) , showline=T  ),
         yaxis = list (title = colnames(data)[2] , gridwidth=2, zeroline=F)
         
    )

}

3/ Determine LOD threshold

We used R-QTL to do some permutation and determine which is the LOD threshold over which q QTL can be considered as significant. (This is done on a cluster because of long computation time). A lod threshold was calculated for each trait. The highest treshold observed was 3.73. We thus used this threshold for every trait since it is the most conservative one.

cd ~/work/FUSA/LOD_SEUIL
cd /NAS/g2pop/HOLTZ_YAN_DATA/DIC2_SILUR/QTL/PUBLI
# je créé un script avec:
for i in $(seq 2 318) ; do echo $i ; Rscript /NAS/g2pop/HOLTZ_YAN_DATA/programmes/Calculate_LOD_threshold.R genotypage.csv phenotypage.csv carte $i ; done
# puis j'envoie
qsub -q normal.q -b yes -cwd -N LOD_thres "./script"
more LOD_thres.o* | grep -B1 -A2 "1000 permutations)" | grep -v "1000 perm" | grep -v "lod" | sed 's/\[1\] \"//' | sed 's/\"//' | sed 's/5% //' > LOD_threshold_per_trait
cat LOD_threshold_per_trait | grep -v "-" | grep -v "[A-Z]" > tmp
R
data=read.table("tmp")
max(data[,1])
min(data[,1])
hist(data[,1])

4/ Agronomical features

Let’s check if we found the expected QTLs for height and flowering date (already known in the litterature..).

height

Supposed to be on chromosome 4B. As you can see on the graphic below, there is no doubt concerning this QTL. It also proove that scripts are working properly. Note you can zoom and hover the graphics to get more details!

plot_a_interactiveQTL(data, "4B", HEI  , my_QTL_thres, 40)

Same result in a table:

summary_QTL_HEI=summary_table_for_QTL(data, HEI, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_HEI))
pop carac chromo LOD_max position marker IC R2 a diff
DS BEQ11_HEI 4B 10.07464 57 Cluster_3278|Contig2|original@188 57-57 11.33570 6.912055 Dic2>Silur
DS CAP13_HEI_blup 4B 21.42926 57 Cluster_3278|Contig2|original@188 57-57 36.53591 0.000000 Dic2>Silur
DS GRI13_HEI_blup 4B 67.63065 57 TRIDC4BG006730.1@987 57-57 64.86658 0.000000 Dic2>Silur
DS GRI15_HEI_blup 4B 11.73699 57 Cluster_3278|Contig2|original@188 57-57 24.19335 0.000000 Silur>Dic2

We can see that this QTL is really consistent from one experiment to another. Dic2 is taller than Silur. Let’s visualize it in a boxplot:

boxplot_a_QTL("TRIDC4BG006730.1@987", "CAP13_HEI_blup"  )

flowering (heading?) date

Flowering date and precocity are known to be controlled by a QTL on chromosome 7B. We also find this result:

p1=plot_a_interactiveQTL(data, "7B", EPI  , my_QTL_thres, 30)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
p2=plot_a_interactiveQTL(data, "7B", FLO  , my_QTL_thres, 30)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
subplot(p1, p2)

Verse

There is only one QTL for the “verse”, on the same position that the QTL that controls height. This is quite logical: the higher is the plant, the more chance it has to fall… Verse is available for CAP13 only. It is funny to see that the QTL is significant for rep1 and rep2, but not for the blup…

plot_a_interactiveQTL(data, "4B", VER , my_QTL_thres, 30)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Epaisseur des glumes

TODO

5/ Resistance to fusariose

Toxine concentration (DON)

Let’s summarize in a table every significant QTL for the “DON” variable:

summary_QTL_DON=summary_table_for_QTL(data, DON, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_DON))
pop carac chromo LOD_max position marker IC R2 a diff
DS CAP13_DON_blup 3A 3.930642 168.5 TRIDC3AG066720.1@1104 164.4-175.3 10.4574 95.60727 Silur>Dic2

If we consider only blups, there is only one QTL on chromosome 3A, with a LOD of 4.19 for DON in CAP13.
Let’s show this 3A chromosome with the DON variables:

plot_a_interactiveQTL(data, "3A", intersect(DON,PUBLI), my_QTL_thres, 6)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

There is absolutely no effect for other trials than CAP13…

Type 1: Note FUSA

Summary of all the significant QTLs for the ** visual notation**:

summary_QTL_NOT=summary_table_for_QTL(data, NOT, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_NOT))
pop carac chromo LOD_max position marker IC R2 a diff
DS BAR14_NOT600_blup 1B 4.976054 1.8 TRIDC1BG001890.1@27 0-13.2 11.5233980 0.1354486 Dic2>Silur
DS BEQ11_NOT450 1B 6.309231 0.5 TRIDC1BG001340.4@2565 0.5-5.1 14.5990816 0.6443723 Dic2>Silur
DS BEQ11_NOT500 1B 7.822870 0.5 TRIDC1BG001340.4@2373 0.5-1.8 17.7999965 0.6567100 Dic2>Silur
DS BEQ11_NOTAUDPC 1B 4.893092 0.5 Cluster_1882|Contig7|original@596 0-6.2 0.9858993 106.7617271 Dic2>Silur
DS GRI11_NOT450_blup 1B 3.936324 142.9 Cluster_16753|Contig2|original@619 136.8-146.4 10.2581917 0.1430546 Silur>Dic2
DS GRI13_NOT550_blup 1B 4.254277 160.4 TRIDC1BG068250.1@1950 160.4-173.9 9.1647534 0.1740732 Silur>Dic2
DS GRI15_NOT350_blup 1B 6.450241 5.1 TRIDC0UG001570.3@1699 1.3-14.5 14.4900190 0.0673185 Dic2>Silur
DS GRI15_NOT550_blup 1B 7.803936 1.3 TRIDC1BG000060.1@264 0.5-6.2 17.3751226 0.2101685 Dic2>Silur
DS GRI15_NOTAUDPC_blup 1B 8.064012 5.1 TRIDC0UG001570.3@2193 1.3-6.2 17.8945425 43.2173606 Dic2>Silur
DS GRI13_NOT350_blup 3A 4.777256 183.7 TRIDC3AG069520.1@1240 175.3-195.5 10.3227644 0.0691291 Silur>Dic2
DS GRI13_NOT450_blup 3A 4.139495 183.7 TRIDC3AG069520.1@672 173.7-190 8.9076580 0.0777861 Silur>Dic2
DS GRI13_NOTAUDPC_blup 3A 5.142387 183.7 TRIDC3AG069520.1@672 175.3-190 11.1175468 53.7728680 Silur>Dic2
DS BEQ11_NOT300 5A 6.468229 239.3 Cluster_265|Contig1|likelySeq@300 236.9-239.3 11.3610513 0.3236671 Silur>Dic2
DS BEQ11_NOT350 5A 7.519043 237.5 Cluster_149|Contig3|original@533 233.9-239.3 17.1751744 0.5422521 Silur>Dic2
DS BEQ11_NOT450 5A 5.981729 247.2 TRIDC5AG069790.7@31 237.5-248.3 13.8769478 0.6393258 Silur>Dic2
DS BEQ11_NOTAUDPC 5A 7.708366 237.5 TRIDC5AG067550.4@2885 236.9-247.2 1.5192131 146.0311422 Silur>Dic2
DS GRI15_NOT550_blup 7B 4.645063 166.2 TRIDC7BG065740.3@1878 157.1-166.2 10.5634901 0.1409060 Dic2>Silur

I have 2 main QTLs: on the 1B and on the 5A.
Then I have some really small messages on 3A for Gri13, and 7B for Gri15

Let’s show the QTLs on graphics. We show only the last note and audpc for blup when available.

toshow=c("GRI11_NOT550_blup","GRI11_NOTAUDPC_blup","GRI13_NOT550_blup","GRI13_NOTAUDPC_blup","GRI15_NOT550_blup","GRI15_NOTAUDPC_blup","BEQ11_NOT500","BEQ11_NOTAUDPC","BAR14_NOT600_blup")
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", intersect(NOT,PUBLI), my_QTL_thres, 10)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", intersect(NOT,PUBLI), my_QTL_thres, 10)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Let’s visualize an example trough a boxplot:

a=boxplot_a_QTL("TRIDC0UG001570.3@2193", "GRI15_NOTAUDPC_blup" )
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
b=boxplot_a_QTL("TRIDC5AG067550.4@2885", "BEQ11_NOTAUDPC")
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
subplot(a,b)

Type 1: % Spike with Fusa

A lot of significant QTLs.

summary_QTL_PEPI=summary_table_for_QTL(data, PEPI, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_PEPI))
pop carac chromo LOD_max position marker IC R2 a diff
DS BAR14_PEPI350_blup 1B 7.383115 5.1 TRIDC0UG001570.3@2193 0.5-5.1 16.798434 7.1975506 Dic2>Silur
DS BAR14_PEPI400_blup 1B 6.446925 5.1 TRIDC0UG001570.3@1806 0.5-14.5 14.813339 6.9793838 Dic2>Silur
DS BAR14_PEPI550_blup 1B 8.231304 11.9 TRIDC1BG003330.11@3771 0.5-13.2 18.524351 10.6648096 Dic2>Silur
DS BAR14_PEPI600_blup 1B 11.514558 1.8 TRIDC1BG001560.1@330 1.3-1.8 24.614519 2.9335351 Dic2>Silur
DS BAR14_PEPIAUDPC_blup 1B 9.760591 1.3 TRIDC1BG001970.2@434 1.3-6.2 21.473043 3897.5655867 Dic2>Silur
DS GRI13_PEPI450_blup 2B 4.436633 14.6 TRIDC2BG002330.3@1130 12.2-18 9.571098 0.2884493 Dic2>Silur
DS GRI13_PEPIAUDPC_blup 2B 3.788355 14.6 TRIDC2BG002330.3@1130 12.2-18 8.114487 506.7832162 Dic2>Silur
DS GRI13_PEPI350_blup 3A 5.065082 183.7 TRIDC3AG069520.1@672 180.1-184.9 10.950229 1.9672612 Silur>Dic2
DS GRI13_PEPI450_blup 3A 3.899620 184.9 Cluster_1251|Contig2|complementarySeq@569 175.3-190 8.366828 0.2872906 Silur>Dic2
DS GRI13_PEPIAUDPC_blup 3A 4.980949 183.7 TRIDC3AG069520.1@1240 179-190 10.767507 610.7678883 Silur>Dic2
DS BEQ11_PEPI350 4A 4.074615 107.6 Cluster_948|Contig2|original@2756 105.9-107.6 8.783936 6.4687500 Dic2>Silur
DS BEQ11_PEPIAUDPC 4A 4.030530 107.6 TRIDC4AG045960.1@749 105.9-107.6 8.613270 3208.9948529 Dic2>Silur
DS BEQ11_PEPI350 4B 4.112312 48.8 TRIDC4BG005390.2@1159 47.6-69.2 8.868261 7.0819805 Silur>Dic2
DS BAR14_PEPI350_blup 5A 5.890239 229.1 TRIDC5AG065190.1@501 229.1-237.5 13.593976 5.9203939 Silur>Dic2
DS BAR14_PEPI400_blup 5A 5.044484 229.1 TRIDC5AG065260.2@750 229.1-238 11.680486 6.7378080 Silur>Dic2
DS BAR14_PEPI600_blup 5A 4.693709 237.5 TRIDC5AG067550.4@2885 233.9-239.3 10.867117 1.8876568 Silur>Dic2
DS BAR14_PEPIAUDPC_blup 5A 6.422906 229.1 TRIDC5AG065260.2@750 229.1-237.5 14.762259 2854.9994215 Silur>Dic2
DS BEQ11_PEPI350 5A 5.310014 237.5 Cluster_343|Contig2|likelySeq@420 236.9-257.6 11.480355 8.7632477 Silur>Dic2
DS BEQ11_PEPIAUDPC 5A 5.086288 237.5 Cluster_149|Contig5|original@1002 233.9-257.6 10.911406 4406.6703131 Silur>Dic2
DS GRI13_PEPI550_blup 5A 3.916479 252.4 TRIDC5AG071590.5@815 246-277.5 8.404146 1.7975920 Silur>Dic2
DS GRI15_PEPI550_blup 5A 4.707260 242.1 Cluster_11094|Contig1|complementarySeq@313 232.1-245.6 10.644227 3.7144168 Silur>Dic2
DS GRI15_PEPIAUDPC_blup 5A 4.639353 242.1 Cluster_13017|Contig1|original@433 231.3-245.6 10.488543 781.0038900 Silur>Dic2
DS CAP13_PEPIAUDPC_blup 6A 3.929697 157.5 TRIDC6AG056520.1@2331 157.5-163.3 10.311553 826.7490335 Silur>Dic2
DS BEQ11_PEPI450 7B 4.395495 39.0 TRIDC7BG009820.4@918 16.8-47.2 9.192045 8.5826325 Dic2>Silur
DS CAP13_PEPI400_blup 7B 8.080811 16.8 TRIDC7BG001540.1@113 16.8-16.8 17.034043 11.3752637 Silur>Dic2
DS GRI15_PEPI550_blup 7B 4.684315 157.9 Cluster_6415|Contig2|original@249 157.1-166.2 10.591689 3.4147478 Dic2>Silur
DS GRI15_PEPIAUDPC_blup 7B 3.975488 166.2 TRIDC7BG065740.3@489 157.1-168.4 8.944878 622.8496833 Dic2>Silur

Two main QTLs on chromosome 1B and 5A. For the 1B QTLs, it is interesting to see that once more, the resistance allele comes from Silur.. How is it possible? The 1B QTLs works for BAR14 only, but for every traits and with strong LOD-scores.The 5A works for 3 experiments.
Then we have some small signals on chromosomes 2A, 2B, 3A, 4A, 6A. Nothing really interesting.
And once more the QTL linked with height and precocity on chromosomes 4B and 7B. For the 4B–> no BLUP.

We can have a look to the LOD-scores of the main traits (AUDPC and LAST):

# Variable to show on the plot?
my_var_to_show=c( Reduce(intersect, list(PEPI,BLUP,AUDPC)), Reduce(intersect, list(PEPI,BLUP,LAST))  )
par(mfrow=c(2,2))
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", my_var_to_show, my_QTL_thres, 12)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", my_var_to_show, my_QTL_thres, 7)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

An example for QTL 1B and 5A on a boxplot:

Let’s visualize an example trough a boxplot:

a=boxplot_a_QTL("TRIDC1BG001970.2@434", "BAR14_PEPIAUDPC_blup")
b=boxplot_a_QTL("TRIDC5AG065260.2@750", "BAR14_PEPIAUDPC_blup")
subplot(a,b)

Type 2: % of Spikelet with Fusa

A lot of significant QTLs.

summary_QTL_PEPIL=summary_table_for_QTL(data, PEPIL, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_PEPIL))
pop carac chromo LOD_max position marker IC R2 a diff
DS BAR14_PEPIL300_blup 1B 5.082209 1.3 TRIDC1BG001970.2@352 0.5-6.2 11.768093 0.9102612 Dic2>Silur
DS BAR14_PEPIL400_blup 1B 5.209982 1.8 TRIDC0UG002680.1@311 0-14.5 12.061141 1.6308273 Dic2>Silur
DS BAR14_PEPIL550_blup 1B 10.983280 1.3 TRIDC1BG000060.1@264 0.5-5.1 23.688770 5.1209416 Dic2>Silur
DS BAR14_PEPIL600_blup 1B 47.428998 1.3 TRIDC1BG001970.2@582 1.3-1.3 58.793684 20.7970099 Dic2>Silur
DS BAR14_PEPILAUDPC_blup 1B 15.252221 1.3 TRIDC1BG001970.2@502 1.3-1.8 30.560263 1499.2994809 Dic2>Silur
DS CAP13_PEPIL300_blup 3B 3.869784 128.6 Cluster_312|Contig3|original@253 108.7-131.6 8.907305 0.9340312 Silur>Dic2
DS BAR14_PEPIL300_blup 5A 3.845809 237.5 Cluster_343|Contig2|likelySeq@420 219.6-247.2 8.850201 0.8289052 Silur>Dic2
DS BAR14_PEPIL550_blup 5A 3.915950 229.1 TRIDC5AG065260.2@750 220.6-246 9.019606 2.8962405 Silur>Dic2
DS CAP13_PEPIL400_blup 7B 5.319279 16.8 TRIDC7BG001540.1@113 16.8-18.7 12.866356 4.8895228 Silur>Dic2

QTL 1B is present for BAR14 only, once more with the Silur allele beeing resistant. The 5A QTL is once more available, with the Dic2 allele beeing resistant. Small QTL on 3B. 4B and 7B QTL are present as usual but nor for BLUPs!

We can have a look to the LOD-scores:

# Variable to show on the plot?
my_var_to_show=c( Reduce(intersect, list(PEPIL,BLUP,AUDPC)), Reduce(intersect, list(PEPIL,BLUP,LAST))  )
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", my_var_to_show, my_QTL_thres, 45)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", my_var_to_show, my_QTL_thres, 7)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Type 2: Nb of Spikelet with Fusa

A lot of significant QTLs.

summary_QTL_NEPIL=summary_table_for_QTL(data, NEPIL, my_QTL_thres, 1.5)
kable(xtable(summary_QTL_NEPIL))
pop carac chromo LOD_max position marker IC R2 a diff
DS BAR14_NEPIL300_blup 1B 4.450629 1.8 TRIDC0UG002680.1@315 0-11.9 10.295939 4.8652348 Dic2>Silur
DS BAR14_NEPIL400_blup 1B 4.751407 1.3 TRIDC1BG001970.2@458 0.5-14.5 11.001903 8.6271445 Dic2>Silur
DS BAR14_NEPIL550_blup 1B 10.625259 1.3 TRIDC1BG000060.1@264 0.5-5.1 23.052405 26.6781814 Dic2>Silur
DS BAR14_NEPIL600_blup 1B 47.399229 1.3 TRIDC1BG000060.1@264 1.3-1.3 58.778490 102.8526378 Dic2>Silur
DS BAR14_NEPILAUDPC_blup 1B 13.706392 1.8 TRIDC0UG002680.1@315 1.3-1.8 28.215512 7778.0452855 Dic2>Silur
DS BEQ11_NEPIL450 1B 9.653514 1.3 TRIDC1BG001970.2@596 0.5-1.8 5.980977 37.0416667 Dic2>Silur
DS BEQ11_NEPIL500 1B 7.761044 1.3 TRIDC1BG001970.2@457 0.5-1.8 2.257181 34.1111111 Dic2>Silur
DS BEQ11_NEPILAUDPC 1B 7.964935 1.8 TRIDC1BG001970.2@423 0.5-5.1 4.195497 4895.1983447 Dic2>Silur
DS GRI13_NEPIL550_blup 1B 3.791953 173.9 TRIDC1BG070800.4@550 160.4-189.4 8.122491 9.3823938 Silur>Dic2
DS GRI11_NEPIL550_blup 2B 4.144849 108.5 Cluster_18588|Contig1|original@301 98.6-110 10.816610 3.7106381 Dic2>Silur
DS BEQ11_NEPIL300 5A 8.214926 237.5 Cluster_149|Contig3|original@533 237.5-239.3 17.186178 6.5306102 Silur>Dic2
DS BEQ11_NEPIL350 5A 7.324755 237.5 Cluster_149|Contig5|original@1002 233.9-239.3 8.609625 18.0060217 Silur>Dic2
DS BEQ11_NEPIL450 5A 4.687992 245.6 Cluster_6432|Contig1|original@938 237.5-250.1 3.053695 28.2070833 Silur>Dic2
DS BEQ11_NEPILAUDPC 5A 6.987653 247.2 TRIDC5AG069790.7@31 237.5-248.3 3.725141 4648.8848315 Silur>Dic2
DS GRI11_NEPIL450_blup 5A 4.005286 46.0 TRIDC5AG005440.2@840 39.7-47.6 10.443455 0.6073008 Silur>Dic2
DS GRI11_NEPIL500_blup 5A 3.854576 47.6 Cluster_7428|Contig3|original@3241 36.4-53.4 10.037839 1.1677623 Silur>Dic2
DS BAR14_NEPIL_blup 7A 5.877454 185.3 Cluster_2718|Contig1|original@181 185.3-208.3 13.563564 12.0307361 Dic2>Silur
DS CAP13_NEPIL400_blup 7B 5.054352 16.8 TRIDC7BG001540.1@113 16.8-18.7 12.235941 8.8491561 Silur>Dic2

The 1B QTL is present for BAR14 as usual, with Silur as resistant. Warning, GRI13 as a QTL on the 1B also, but not on the same position.
Small QTL on 1A, 2B, 3A, 3B, 4A, 6B, 7A. QTL height and preco as usual, but not for BLUP.

We can have a look to the LOD-scores:

# Variable to show on the plot?
my_var_to_show=c( Reduce(intersect, list(NEPIL,BLUP,AUDPC)), Reduce(intersect, list(NEPIL,BLUP,LAST))  )
par(mfrow=c(2,2))
#Chromosome 1B
plot_a_interactiveQTL(data, "1B", my_var_to_show, my_QTL_thres, 45)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
#Chromosome 5A
plot_a_interactiveQTL(data, "5A", my_var_to_show, my_QTL_thres, 7)
## Warning: We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

6/ Summary table

Per experiment

Let’s summarize what QTL has been found for each experiment. For each expe and each phenotype, we write QTL found, and the number of phenotypes that are significant for this QTL between bracklets.

# initialize table
bil=data.frame(matrix(0,5,7))
colnames(bil)=c("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11", "BAR14")
rownames(bil)=c("DON", "NOT","PEPI", "NEPIL", "PEPIL")
# fill it
# For each variable (DON, EPIL,...)
num=0
for(j in list(summary_QTL_DON, summary_QTL_NOT, summary_QTL_PEPI, summary_QTL_NEPIL, summary_QTL_PEPIL)){
  num=num+1
  
  # And for each experiment
  num2=0
  for(i in colnames(bil)){
     num2=num2+1
     a=j[ grep(i, j$carac ) , ]
     b=table(a$chromo)
     c=paste(names(b), " (", b, ")", sep="")
     d=paste(c, collapse = " - ")
     bil[num,num2]=d
  }}

# write when data is not available
bil[1,c(2,4,6,7)]="no data"
bil[2,c(1)]="no data"
bil[3,c(2)]="no data"
bil[5,c(2,6)]="no data"
CAP13 GRI11 GRI13 GRI15 LEC14 BEQ11 BAR14
DON 3A (1) no data () no data () no data no data
NOT no data 1B (1) 1B (1) - 3A (3) 1B (3) - 7B (1) () 1B (3) - 5A (4) 1B (1)
PEPI 6A (1) - 7B (1) no data 2B (2) - 3A (3) - 5A (1) 5A (2) - 7B (2) () 4A (2) - 4B (1) - 5A (2) - 7B (1) 1B (5) - 5A (4)
NEPIL 7B (1) 2B (1) - 5A (2) 1B (1) () () 1B (3) - 5A (4) 1B (5) - 7A (1)
PEPIL 3B (1) - 7B (1) no data () () () no data 1B (5) - 5A (2)

Details of all QTL

Let’s save a huge table that will contain absolutely every QTL found as significant in this study. It will be a supplementary data for the publication.

OR=rbind(summary_QTL_DON, summary_QTL_NOT, summary_QTL_PEPI, summary_QTL_NEPIL, summary_QTL_PEPIL)

All signif QTL for publication

In the publication we will talk about QTL on blups only, when available, and on raw data if not. Let’s list of these QTLs:

OR=OR[which(OR$carac%in%PUBLI) , ]
dim(OR)
## [1] 72 10
OR=OR[order(OR$chromo, OR$carac), ]
dim(OR)
## [1] 72 10
OR$position=round(OR$position,2)
write.table(OR, file="~/Dropbox/Publi_Fusariose/SUPPORTING_DATA/OR_detail_all_QTLs.csv", row.names=F, col.names=T, quote=F, sep=";", dec=,)

Save environment

save(list=ls() , file="/Users/yan/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/QTL_R_environment.R")

Yan Holtz

October 2016